library("IRdisplay")
# Enseño el grafo dado por el experto
display_png(file="grafo.png")
NOTA: Como ya se comentado anteriormente, es importante evitar tildes a la hora de nombrar los nodos. Tampoco son recomendables espacios en blanco ni ningún tipo de carácter especial.
A continuación, introduce el DAG en R utilizando la definicíon del modelo que acabas de crear
Realiza una lista de los padres e hijos de cada uno de los nodos
Realiza una lista de todas las conexiones fundamentales presentes en el grafo, y realiza una clasificación de cada una de ellas determinando si es una estructura en serie, divergente o convergente.
¿Se observa alguna v-estructura en el grafo?
Intenta introducir un arco que parta del nodo neblina y se dirija hacia escarcha, y otro que vaya de granizo a lluvia. Comenta qué sucede en cada caso, y si el resultado daría lugar a una red bayesiana válida.
Determina la manta de Markov del nodo rocio
Introduce un arco que parta del nodo Lluvia y se dirija al nodo Niebla, ¿cambia la manta de Markov del nodo rocio? En caso afirmativo, ¿cómo lo hace?
Para continuar, elimina el arco definido anteriormente entre los nodos Lluvia y Niebla, recuperando el DAG original.
[Viento][Rocio][Escarcha][Niebla][Nieblina][Lluvia][Granizo][Tormenta][Nieve][Nieve_Suelo]
library(bnlearn)
## Defino el grafo vacío:
dag <- empty.graph(nodes=c("viento", "rocio", "escarcha", "niebla", "neblina",
"lluvia", "granizo", "tormenta", "nieve", "nieveSuelo"))
# Defino los arcos
arc.set <- matrix(c("viento", "lluvia",
"viento", "niebla",
"viento", "rocio",
"rocio", "niebla",
"rocio", "escarcha",
"niebla", "neblina",
"lluvia", "tormenta",
"tormenta", "granizo",
"tormenta", "nieve",
"nieve", "nieveSuelo"),
byrow = TRUE, ncol = 2,
dimnames = list(NULL, c("from", "to")))
arcs(dag) <- arc.set
print(dag)
plot(dag)
# Verifico la factorización del punto anterior
modelstring(dag)
for (node in nodes(dag)){
print(paste("NODO: ", node))
print("----")
print("PADRES:")
print(parents(dag, node = node))
print("----")
print("HIJOS:")
print(children(dag, node = node))
print("====")
}
Se encuentra la estructura convergente
viento --> niebla <-- rocio
Por otro lado, esa no representa una v-estrucura, porque hay un arco que conecta viento con rocio.
El grafo presenta 3 estructuras divergentes.
Estructuras en serie:
No hay v-estructuras en el grafo
# Verifico que no haya v-estructuras
vstructs(dag)
dag <- set.arc(dag, from = "neblina", to = "escarcha")
vstructs(dag)
plot(dag)
dag <- set.arc(dag, from = "granizo", to = "lluvia")
mb(dag, "rocio")
La introducción de este arco genera una nueva v-estrucura: rocio --> niebla <-- lluvia.
Esto hace que el nodo lluvia entre en la manta de Markov del nodo rocio.
dag <- set.arc(dag, from = "lluvia", to = "niebla")
vstructs(dag)
mb(dag, "rocio")
plot(dag)
# Vuelvo a definir los arcos originales
arc.set <- matrix(c("viento", "lluvia",
"viento", "niebla",
"viento", "rocio",
"rocio", "niebla",
"rocio", "escarcha",
"niebla", "neblina",
"lluvia", "tormenta",
"tormenta", "granizo",
"tormenta", "nieve",
"nieve", "nieveSuelo"),
byrow = TRUE, ncol = 2,
dimnames = list(NULL, c("from", "to")))
arcs(dag) <- arc.set
print(dag)
plot(dag)
Considerando la información proporcionada por el conjunto de datos meteor, construye la red bayesiana utilizando el método bayesiano de estimación paramétrica.
¿Cuál sería el número potencial de parámetros (sin usar la red bayesiana) del modelo para calcular la probabilidad global si no utilizasemos el DAG?
¿Cuántos parámetros tiene la distribución global dada por la red bayesiana?
Determina el número de parámetros asociado a cada una de las distribuciones locales
Obtén las tablas de probabilidad condicional asociadas los nodos granizo y niebla. Ahora representa la información de cada tabla en sendos gráficos.
# Construyo la red bayesiana a partir de la tabla
meteoro <- read.table("meteoro.txt", header = TRUE)
bn = bn.fit(dag, data = meteoro, method = "bayes")
número de parametros = p = $2^n$ - 1, con n = numero de nodos = 10
p = $2^{10} - 1$ = 1023
Considerando el grafo, que hemos contruido, los parametros son dados por la factorización:
[viento][rocio|viento][lluvia|viento][escarcha|rocio][niebla|viento:rocio] [tormenta|lluvia][neblina|niebla][granizo|tormenta][nieve|tormenta] [nieveSuelo|nieve]
Miramos cuantos parametros tiene cada uno de los factores:
[viento]: 1 parametro (probabilidad de 's' y probabilidad de 'n' = 1 - p('s'))
[rocio|viento]: 2 parametros
[lluvia|viento]: 2 parametros
[escarcha|rocio]: 2 parametros
[niebla|viento:rocio] : 4 parametros
[tormenta|lluvia]: 2 parametros
[neblina|niebla]: 2 parametros
[granizo|tormenta]: 2 parametros
[nieve|tormenta]: 2 parametros
[nieveSuelo|nieve]: 2 parametros
Sumando: 21
# Verifico con el código
nparams(bn)
# Nodo granizo
print(bn$granizo$prob)
plot(bn$granizo$prob)
# Nodo niebla
print(bn$niebla$prob)
plot(bn$niebla$prob)
Una vez construída la red hemos establecido la base de conocimiento del sistema inteligente. A continuación se puede calcular la probabilidad de cualquier variable o conjunto de variables condicionadas a cualquier evidencia que se tenga disponible para un problema dado, es decir, realizar la inferencia.
NOTA: En este apartado se deberá aplicar la inferencia exacta. Conocido que en un día dado se han producido tormentas, calcula cómo afecta este hecho a la probabilidad de que se produzcan los siguientes fenómenos meteorológicos:
Dado el ejercicio anterior, repetir ahora el ejercicio mediante inferencia aproximada, calculando para cada una de las estimaciones 100 realizaciones. Para cada una de las respuestas anteriores, representa un diagrama de cajas que represente la dispersión en la estimación de la probabilidad, marcando además el valor obtenido de manera exacta en el apartado anterior.
# Vuelvo a pintar el grafo
plot(dag)
La nieve y el granizo son fenómenos independientes a priori
No, ya que tienen el mismo padre y no son d-separados. Podrían ser independientes, pero no lo podemos afirmar con lo que conocemos.
# Verifico
dsep(dag, x = "nieve", y = "granizo")
La nieve y el granizo son fenómenos independientes dado que haya habido tormenta
Sí, en este caso el conocer el estado del nodo tormenta hace que el camino entre nieve y tormenta no esté activo. Por lo tanto, los dos nodos son d-separados y entonces independientes.
# Verifico
dsep(dag, x = "nieve", y = "granizo", z = "tormenta")
La nieve en el suelo y la neblina son fenómenos independientes
No: como existen dos caminos que las conectan, podrían ser independientes, pero no lo podemos afirmar con lo que conocemos.
# Verifico
dsep(dag, x = "nieveSuelo", y = "neblina")
La nieve en el suelo y la neblina son fenómenos independientes dado que haya habido tormenta
Sí, como en caso de nieve y granizo, si conocemos el estado del nodo tormenta se interrumpe el camino entre niebla y nieve en el suelo, y las dos variables se d-separan.
# Verifico
dsep(dag, x = "nieveSuelo", y = "neblina", z = "tormenta")
library(gRain)
# Compilo la red
red <- compile(as.grain(bn))
# Fijo el estado del nodo tormenta a 's'
red.tormenta.s = setEvidence(red, nodes = "tormenta", states = "s")
# Consulto la red: P(lluvia=s|tormenta=s) = ?
query1 <- querygrain(red.tormenta.s, nodes = "lluvia", type = "marginal")$lluvia["s"]
print(paste("La probabilidad de lluvia dado que hay tormenta es", query1))
# Consulto la red: P(viento=s|tormenta=s) = ?
query2 <- querygrain(red.tormenta.s, nodes = "viento", type = "marginal")$viento["s"]
print(paste("La probabilidad de viento dado que hay tormenta es", query2))
# Consulto la red: P(lluvia=s,viento=n|tormenta=s) = ?
query3 <- querygrain(red.tormenta.s, nodes = c("lluvia","viento"), type = "joint")["s","s"]
print(paste("La probabilidad de viento dado que hay tormenta es", query3))
# Miro la probabilidad de que llueva, sin poner condiciones
query4 <- querygrain(red, nodes = "lluvia", type = "marginal")$lluvia["s"]
print(paste("La probabilidad de lluvia dado que hay tormenta es", query1,
"mientras que sin poner condiciones, la probabilidad de que llueva es", query4,
". Eso quiere decir que es más probable que llueva cuando hay tormenta."))
# Miro la probabilidad de que haya viento, sin poner condiciones
query5 <- querygrain(red, nodes = "viento", type = "marginal")$viento["s"]
print(paste("La probabilidad de que haya viento dado que hay tormenta es", query2,
"mientras que sin poner condiciones, la probabilidad de que haya viento es", query5,
". Eso quiere decir que es más probable que haya viento cuando hay tormenta.",
"Cuantitativamente, el aumento es del", 100*(query2 / query5 - 1), "%"))
# Fijo el numero de simulaciones a 100
N = 100
# Simulo 100 veces lo que pasa a la red sin poner condiciones
bn.sim <- simulate(red, n = N)
# Verifico que tormenta tenga como valores tantos 's' como 'n'
print(paste("Número de 's'", nrow(bn.sim[bn.sim$tormenta == 's',])))
print(paste("Número de 'n'", nrow(bn.sim[bn.sim$tormenta == 'n',])))
# Miro la probabilidad de que haya tormenta, sin poner condiciones
# para verificar que los números que saco tienen sentido
query.tormenta <- querygrain(red, nodes = "tormenta", type = "marginal")$tormenta
query.tormenta
# Simulo 100 veces lo que pasa a la red si siempre hay tormenta
bn.sim.tormenta.s <- simulate(red.tormenta.s, n = 100)
# Verifico que tormenta siempre tenga como valor 's'
bn.sim.tormenta.s$tormenta
Ahora que tengo las simulaciones, puedo hacer consultas:
prob.lluvia <- nrow(bn.sim.tormenta.s[bn.sim.tormenta.s$lluvia == 's',])/N
print(paste("Probabilidad de lluvia, dada tormenta:",
prob.lluvia))
my.table = table(bn.sim.tormenta.s$lluvia)
plot(my.table, type = "p", ylim = c(0,100))
arrows(1, my.table["n"] + sqrt(my.table["n"]),
1, my.table["n"] - sqrt(my.table["n"]),
length=0.05, angle=90, code=3)
points(1, 100*(1 - query1), col = "red")
arrows(2, my.table["s"] + sqrt(my.table["s"]),
2, my.table["s"] - sqrt(my.table["s"]),
length=0.05, angle=90, code=3)
points(2, 100*query1, col = "red")
prob.viento <- nrow(bn.sim.tormenta.s[bn.sim.tormenta.s$viento == 's',])/N
print(paste("Probabilidad de viento, dada tormenta:",
prob.viento))
my.table = table(bn.sim.tormenta.s$viento)
plot(my.table, type = "p", ylim = c(0,100))
arrows(1, my.table["n"] + sqrt(my.table["n"]),
1, my.table["n"] - sqrt(my.table["n"]),
length=0.05, angle=90, code=3)
points(1, 100*(1 - query2), col = "red")
arrows(2, my.table["s"] + sqrt(my.table["s"]),
2, my.table["s"] - sqrt(my.table["s"]),
length=0.05, angle=90, code=3)
points(2, 100*query2, col = "red")
Que llueva y que además las rachas de viento superen los 50 Km/h, P(lluvia=s,viento=s|tormenta=s)
prob.viento.lluvia <- nrow(bn.sim.tormenta.s[bn.sim.tormenta.s$viento == 's' & bn.sim.tormenta.s$lluvia == 's',])/N
print(paste("Probabilidad de viento y lluvia, dada tormenta:",
prob.viento.lluvia))
my.table = table(bn.sim.tormenta.s$viento == 's' & bn.sim.tormenta.s$lluvia == 's')
names(my.table) = c("n", "s")
plot(my.table, type = "p", ylim = c(0,100))
arrows(1, my.table["n"] + sqrt(my.table["n"]),
1, my.table["n"] - sqrt(my.table["n"]),
length=0.05, angle=90, code=3)
points(1, 100*(1 - query3), col = "red")
arrows(2, my.table["s"] + sqrt(my.table["s"]),
2, my.table["s"] - sqrt(my.table["s"]),
length=0.05, angle=90, code=3)
points(2, 100*query3, col = "red")
print(paste("La probabilidad de lluvia dado que hay tormenta es",
prob.lluvia))
prob.lluvia.no_cond <- nrow(bn.sim[bn.sim$lluvia == 's',])/N
print(paste("Mientras que sin poner condiciones, la probabilidad de que llueva es",
prob.lluvia.no_cond))
print("Eso quiere decir que es más probable que llueva cuando hay tormenta.")
my.table = table(bn.sim$lluvia)
plot(my.table, type = "p", ylim = c(0,100))
arrows(1, my.table["n"] + sqrt(my.table["n"]),
1, my.table["n"] - sqrt(my.table["n"]),
length=0.05, angle=90, code=3)
points(1, 100*(1 - query4), col = "red")
arrows(2, my.table["s"] + sqrt(my.table["s"]),
2, my.table["s"] - sqrt(my.table["s"]),
length=0.05, angle=90, code=3)
points(2, 100*query4, col = "red")
print(paste("La probabilidad de que haya viento dado que hay tormenta es",
prob.viento))
prob.viento.no_cond <- nrow(bn.sim[bn.sim$viento == 's',])/N
print(paste("Mientras que sin poner condiciones, la probabilidad de que haya viento es",
prob.viento.no_cond))
print("Eso quiere decir que es menos probable que haya viento cuando hay tormenta.")
print(paste("Cuantitativamente, disminuye del",
100 * (prob.viento / prob.viento.no_cond - 1), "%"))
my.table = table(bn.sim$viento)
plot(my.table, type = "p", ylim = c(0,100))
arrows(1, my.table["n"] + sqrt(my.table["n"]),
1, my.table["n"] - sqrt(my.table["n"]),
length=0.05, angle=90, code=3)
points(1, 100*(1 - query5), col = "red")
arrows(2, my.table["s"] + sqrt(my.table["s"]),
2, my.table["s"] - sqrt(my.table["s"]),
length=0.05, angle=90, code=3)
points(2, 100*query5, col = "red")
Como hemos visto, es posible realizar un aprendizaje automático de la estructura del DAG a partir de los datos, usando algoritmos específicos para ello. Durante las clases hemos visto el ejemplo del algoritmo hill-climbing, aunque como vimos, hay otras posibilidades. Tambien hemos visto que podemos combinar nuestra experiencia y el aprendizaje automático definiendo previamente relaciones entre variables que queremos introducir o descartar en el DAG resultante. Además, se ha explicado que existen scores que sirven como criterio para evaluar la fuerza de la dependendencia entre nodos de la red y comparar la bondad de ajuste del modelo.
Evalúa la significación de los arcos dibujados por el experto en el actual DAG utilizando el estadístico χ2 ¿Hay algún arco no significativo? ¿Cuáles son los tres pares de nodos que presentan un arco de unión más fuerte?
Además del algoritmo hill-climbing, existe otro algoritmo popular de tipo “voraz” denominado Tabu search. En bnlearn se encuentra implementado a traves de la función tabu, y los argumentos de entrada son similares a los vistos para hill-climbing.
Ahora vuelve a aprender de forma automática el DAG usando tabu y hill-climbing, pero imponiendo las siguientes restricciones: 1. Los arcos viento --> lluvia, tormenta --> granizo y nieve --> nieveSuelo deben quedar reflejados en el DAG. 2. Ningún arco debe unir directamente la neblina con el granizo ni la niebla con la tormenta.
# Uso arc.strength para valorar la significación de los arcos
arc.strength(dag, data = meteoro, criterion = "x2")
Se observa que los p-valores de todos los arcos son muy pequeños ($< 10^{-5}$) y por eso podemos considerar que todos son significativos. Dicho de otra manera, los arcos sugeridos por el experto tienen en cuenta relaciones probabilisticas que de verdad existen en el conjunto de datos.
En particular, los arcos más significativo (los que presentan p-valor más pequeño) son:
# Construyo el DAG usando el algoritmo Tabu search
dag.tabu <- tabu(x = meteoro)
dag.tabu
plot(dag.tabu)
# Construyo el DAG usando el algoritmo hill-climbing
dag.hc <- hc(x = meteoro)
dag.hc
plot(dag.hc)
# BIC score para el DAG original
dag.score <- bnlearn::score(dag, data = meteoro)
dag.score
# BIC score para el DAG generado con el algoritmo Tabu search
dag.tabu.score <- bnlearn::score(dag.tabu, data = meteoro)
dag.tabu.score
# BIC score para el DAG generado con el algoritmo hill-climbing
dag.hc.score <- bnlearn::score(dag.hc, data = meteoro)
dag.hc.score
El mejor BIC (más cerca de 0, o más pequeño en valor absoluto) se obtiene con el algoritmo Tabu search.
# graphviz.plot me dá el siguiente problema en local:
graphviz.plot(dag)
# Por esta razón, he dibujado los grafos en el datasciencehub (donde no logro instalar gRain)
# y los pego aquí abajo
print("DAG experto")
display_png(file="DAG.png")
print("DAG tabu")
display_png(file="DAG_Tabu.png")
# Uso arc.strength para valorar la significación de los arcos
arc.strength(dag.tabu, data = meteoro, criterion = "x2")
print("DAG hill-climbing")
display_png(file="DAG_hc.png")
# Uso arc.strength para valorar la significación de los arcos
arc.strength(dag.hc, data = meteoro, criterion = "x2")
Sin duda el grafo más legible e interpretable es el que ha proporcionado el experto, donde se notan relaciones que podemos considerar 'lógicas' y facilmente explicables, como:
Por otro lado, este conocimiento, aunque importante y fundamental, puede resultar limitado a la sola experiencia y capacidad de observación de quién estudia el problema, mientras que relaciones importantes entre las caracteristicas de los datos pueden estar escondidas y aparecer solo gracias al utilizo de algoritmos de aprendizaje.
Se observa por ejemplo, que en los dos grafos construidos con algoritmos de aprendizaje, aparece un arco (en un caso muy fuerte) entre nieve y granizo, que no aparece en el grafo experto.
Hay que considerar que la otra cara de este tipo de aprendizaje basado solo en los datos puede considerar relaciones probabilisticas debiles y probablemente basadas más en fluctuaciones estadisticas que en reales dependencia entre dos eventos. Por ejemplo, me resulta extraño que haya relación entre granizo rocio. Se pueden generar así arcos que en vez de describir relaciones entre las caracteristicas del problema, los hacen sobre caracteristicas del dataset (como un sobreajuste), haciendo que el número de parametros del modelo crezca sin aportar realmente nada a su capacidad predictiva. Para obviar a este inconveniente, se suelen introducir términos de regularización en la computación del score del modelo, para penalizar modelos demasiado complejos.
# Preparo la whitelist
whitelist <- matrix(c("viento", "lluvia",
"tormenta", "granizo",
"nieve", "nieveSuelo"),
ncol = 2,
dimnames = list(NULL, c("from", "to")))
# Preparo la blacklst
blacklist <- matrix(c("neblina", "granizo",
"granizo", "neblina",
"tormenta", "niebla",
"niebla", "tormenta"),
ncol = 2, byrow = TRUE,
dimnames = list(NULL, c("from", "to")))
# Construyo el DAG "semiexperto" usando el algoritmo Tabu search
dag.tabu.semiexperto <- tabu(x = meteoro,
whitelist = whitelist,
blacklist = blacklist)
dag.tabu.semiexperto
plot(dag.tabu.semiexperto)
dag.tabu.score.semiexperto <- bnlearn::score(dag.tabu.semiexperto, data = meteoro)
dag.tabu.score.semiexperto
# Uso arc.strength para valorar la significación de los arcos
arc.strength(dag.tabu.semiexperto, data = meteoro, criterion = "x2")
# Construyo el DAG "semiexperto" usando el algoritmo hill-climbing search
dag.hc.semiexperto <- hc(x = meteoro,
whitelist = whitelist,
blacklist = blacklist)
dag.hc.semiexperto
plot(dag.hc.semiexperto)
dag.hc.score.semiexperto <- bnlearn::score(dag.hc.semiexperto, data = meteoro)
dag.hc.score.semiexperto
# Uso arc.strength para valorar la significación de los arcos
arc.strength(dag.hc.semiexperto, data = meteoro, criterion = "x2")
# Preparo una tabla de resumen con todos los resultados de los scores
summary.table <- matrix(c("DAG experto", dag.score,
"DAG Tabu", dag.tabu.score,
"DAG hill-climbing", dag.hc.score,
"DAG Tabu semi-experto", dag.tabu.score.semiexperto,
"DAG hill-climbing semi-experto", dag.hc.score.semiexperto),
ncol = 2, byrow = TRUE,
dimnames = list(NULL, c("DAG", "SCORE")))
summary.table
Como hemos visto anteriormente, el uso de algoritmos de aprendizaje ha permitido mejorar el score BIC que se había obtenido usando solo el conocimiento de un experto.
Vemos ahora que combinar las dos estrategias reduce de manera casi imperceptible la capacidad del grafo de adaptarse a los datos, pero permite evitar que relaciones probabilsticas poco significativas entren en el DAG.
Las relaciones incluidas en la whitelist:
Ya aparecían en los grafos construidos a través de algoritmos de aprendizaje, que pero introducian también relaciones que hemos puesto en la blacklist, compuesta por:
Como podemos ver en las tablas más abajo, la relación probabilistica entre granizo y neblina aparece, aunque poco significativa, mejorando así de manera artificial el score de los grafos, al precio de introducir mayor complejidad y cantidad de parametros. Lo mismo vale para arcos como "granizo --> rocio", que tienen una fuerza muy baja y que podemos entender como poco significativos, y que pero encontramos en los grafos construidos sin ningún tipo de conocimiento experto.
En definitiva, el método que me ha parecido más interesante ha sido el que mezcla algoritmos de aprendizaje y conocimiento previo del problema, por la siguientes razones:
arc.strength(dag.hc, data = meteoro, criterion = "x2")
arc.strength(dag.tabu, data = meteoro, criterion = "x2")